Análisis y pronóstico de la información accidentalidad vehicular en el municipio de Medellín para los años 2014 a 2018

2.3.1. Limpieza y resumen de los datos

De manera inicial se leen los archivos.

Nota al ejecutador de código: para algunos sistemas operativos o versiones de paquetes la limpieza no funciona correctamente, por lo que se recomienda volver a cargar el archivo a partir aquí.

## lectura de archivos

# lista archivos
lista <- list.files(pattern = "^Acc.*.csv", include.dirs = T, recursive = T)

# leer todos los archivos
lista_df <- map(lista, fread, sep = ",", encoding = "UTF-8", colClasses = "c")

# agregar archivos del 2014 a 2018
acc <- bind_rows(lista_df)

# ver cabecera del archivo
head(acc)

Se identifica que se deben hacer las siguientes correcciones en la tabla acc:

• Organizar la columna FECHA en formato ISO 8601..

• Debido a la naturaleza del trabajo, se eliminan los nulos de la columna CLASE y se unifican los tipos de accidente, limpiando tildes y convirtiendo en mayúscula.

• De forma similar, se eliminan los nulos de la columna DISENO, se unifican los tipos de diseño, se limpian tildes y se convierte a mayúscula.

• Se eliminan las tildes de los días de la semana en la columna DIA_NOMBRE.

• Se crea la columna COMUNA_BARRIO a partir de la columna CBML con el objetivo de que sirva como clave foránea para la unión con el archivo espacial de barrios, disponible aquí.

## limpieza de la tabla

# organizar fecha
acc$FECHA <- ymd(gsub(pattern = "T.*", replacement = "", acc$FECHA),
                 "%Y-%M-%D")[1:209426]

# eliminar datos nulos y corregir clase
acc <- acc[-which(acc$CLASE == ""),]
acc$CLASE <- iconv(acc$CLASE, from = "UTF-8", to = "ASCII//TRANSLIT")
acc$CLASE <- gsub("DE ", "", toupper(acc$CLASE))

# eliminar datos nulos y corregir disenio
acc <- acc[-which(acc$DISENO == ""),]
acc$DISENO <- iconv(acc$DISENO, from = "UTF-8", to = "ASCII//TRANSLIT")
acc$DISENO <- gsub("DE ", "", toupper(acc$DISENO))

# corregir tildes de DIA_NOMBRE
acc$DIA_NOMBRE <- iconv(acc$DIA_NOMBRE, from = "UTF-8", to = "ASCII//TRANSLIT")

# crear columna de comuna_barrio
acc <- mutate(acc, COMUNA_BARRIO = str_sub(CBML, 1, 4))

# visualizar nueva acc
head(acc)

Se procede a cargar el archivo espacial de barrios para identificar correctamente el barrio y la comuna de la ubicación del accidente. Luego, se realiza una nueva limpieza de la tabla acc:

• Se limpian las tildes de la columna barrio y se convierte a mayúscula.

• Se eliminan los registros duplicados generados por la unión.

• Se eliminan algunas columnas no necesarias para el análisis.

• Se renombran las columnas de NOMBRE_BAR y NOMBRE_COM por BARRIO y COMUNA respectivamente.

## carga de archivo espacial de barrios

# cargar archivo shp de barrios de medellin
barrio <- read_sf("files/Limite_Barrio_Vereda_Catastral/Limite_Barrio_Vereda_Catastral.shp")

# unir columnas de nombre barrio y comuna
acc <- inner_join(acc, 
                  select(barrio, CODIGO, NOMBRE_COM, NOMBRE_BAR),
                  by = c("COMUNA_BARRIO" = "CODIGO"))
  
# limpiar nombre de barrios
acc$NOMBRE_BAR <- iconv(acc$NOMBRE_BAR, from = "UTF-8", to = "ASCII//TRANSLIT")
acc$NOMBRE_BAR <- toupper(acc$NOMBRE_BAR)

# eliminar posibles duplicados por errores de union
acc <- data.table:::unique.data.table(acc, by = "RADICADO")

# eliminar columnas
acc <- select(acc, -BARRIO, -COMUNA, -OBJECTID, -RADICADO, -DIRECCION_ENC,
              -DIRECCION, -HORA, -CBML, -TIPO_GEOCOD, -MES_NOMBRE, -geometry)

# renombrar
names(acc)[12:13] <- c("COMUNA","BARRIO")

# visualizar acc
head(acc)

Se realiza la lectura de la tabla de días especiales, se limpia y se realiza unión con la tabla de accidentes. Además, se cambia las columnas tal que si es un día especial tenga la marcación Si, de lo contrario será No.

## carga de archivo de fechas especiales

# lectura del archivo
festivos <- fread("files/festivos_y_especiales.csv", header = T)

# convertir a si o no
festivos[,names(festivos)[-1]] <- festivos %>%
                                  transmute_at(c(names(festivos)[-1]),
                                               funs(ifelse(. == "X",
                                                           "Si", 
                                                           ifelse(. == "",
                                                                  "No",
                                                                  .)
                                                           )
                                                    )
                                               )

# convetir a formato fecha
festivos$FECHA <- ymd(festivos$FECHA)[1:172]

# visualizar festivos
head(festivos)

Se realiza la unión de las tablas y se convierten los datos de los días especiales a No en caso de que no haya encontrado coincidencia en la unión.

# union las tablas de dias especiales y accidentes

# union de tablas
acc <- merge(x = acc, y = festivos, by = "FECHA", all.x = T)

# transformar variables
acc[,names(festivos)[-1]] <- acc %>%
                             transmute_at(c(names(festivos)[-1]),
                                         funs(ifelse(is.na(.),
                                                     "No",
                                                     .)
                                              )
                                         )

# visualizar acc
head(acc)

Se realiza la conversión a la codificación UTF-8 y se guarda el nuevo archivo.

## guardar archivo definitivo 

# convertir a utf8
acc[,2:length(acc)] <- map(.x = acc[,2:length(acc)], .f = enc2utf8)

 # guardar archivo
fwrite(acc, "files/accidentalidad_georreferenciada_completa.csv", sep = ",")

# eliminar archivo y limpiar memoria
rm(acc)
gc(reset = T)

##           used (Mb) gc trigger  (Mb) max used (Mb)
## Ncells 1457577 77.9    2695477 144.0  1457577 77.9
## Vcells 8397131 64.1   31568764 240.9  8397131 64.1

## carga de archivo definitivo

# volver a cargar el archivo
acc <- fread("files/accidentalidad_georreferenciada_completa.csv",
             sep = ",",
             encoding = "UTF-8")

# organizar fecha
acc$FECHA <- ymd(acc$FECHA)

# visualizar
head(acc)

2.3.2. Análisis

Se realiza un análisis descriptivo para identificar las variables asociadas a la accidentalidad vial en la ciudad de Medellín a través del análisis y la interpretación de las cifras de accidentes registrados en la Secretaría de Movilidad.

Se crea una función para crear gráficos compuestos de varios subgráficos, la cual se llamará multiplot.

## funcion multiplot

#librerias para graficos y subgraficos
library(ggplot2)
library(grid)

# crear funcion multiplot
multiplot <- function(..., plotlist = NULL, file, cols = 1, layout = NULL) {

  # hacer una lista de los argumentos en ... y plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # si layout es NULL, entonces use columnas para establecerlo
  if (is.null(layout)) {
    # crear el panel
    # ncol: numero de columnas del grafico
    # nrow: numero de filas, calculado del numero de columnas
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, 
                     nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # configurar la pagina
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # poner cada grafico en la posicion correcta
    for (i in 1:numPlots) {
      # obtener la posicion i,j de la matriz de la region que va a contener 
      # el subgrafico
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

## grafica dia, mes y anio

# agrupar por nombre del dia
acc_nombredia <- acc %>% 
                 group_by(DIA_NOMBRE) %>% 
                 summarize(total_registros = n())

# establecer orden de los dias
acc_nombredia$DIA_NOMBRE <- ordered(acc_nombredia$DIA_NOMBRE, 
                                    levels = c( "LUNES", "MARTES", "MIERCOLES",
                                                "JUEVES", "VIERNES","SABADO",
                                                "DOMINGO")
                                    )
# grafica dia
p2 <- ggplot(data = acc_nombredia, aes(x = DIA_NOMBRE, y = total_registros)) + 
      geom_bar(stat = "identity", 
               position="dodge", 
               fill = "dodgerblue3", 
               color = "grey48", 
               alpha = 0.8) +
      xlab("Días") + # eje x
      ylab("Total registros") + # eje y
      # Girar label eje x
      theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5)) + 
      ggtitle("Número de accidentes por días de la semana") # titulo del grafico

# agrupar por mes
acc_mes <- acc %>% 
           group_by(MES) %>% 
           summarize(total_registros = n())

# grafica mes
p1 <- ggplot(data = acc_mes, aes(x = MES, y = total_registros)) + 
      geom_bar(stat = "identity", 
               position="dodge", 
               fill = "dodgerblue3", 
               color = "grey48", 
               alpha = 0.8) +
      xlab("Mes") + # eje x
      ylab("Total registros") + # eje y
      ggtitle("Número de accidentes por mes") # titulo del grafico

# agrupar por anio
acc_ano <- acc %>% 
           group_by(PERIODO) %>% 
           summarize(total_registros = n())

# grafica anio
p3 <- ggplot(data = acc_ano, aes(x = PERIODO, y = total_registros)) + 
      geom_bar(stat = "identity", 
               position = "dodge", 
               fill = "dodgerblue3", 
               color = "grey48", 
               alpha = 0.8) +
      xlab("Año") + # eje x
      ylab("Total registros") + # eje y
      ggtitle("Número de accidentes por año") # titulo del grafico

# crear multigrafico
multiplot(p1, p2, p3, cols = 2)

## grafica comuna y disenio

# agrupar por comuna
acc_comuna <- acc %>% 
              group_by(COMUNA) %>% 
              summarize(total_registros = n())

# crear grafica comuna
p4 <- ggplot(data = acc_comuna, aes(x = COMUNA, y = total_registros)) +
      geom_bar(stat = "identity", 
               position = "dodge", 
               fill = "dodgerblue3", 
               color = "grey48", 
               alpha = 0.8) +
      xlab("Comuna") + # eje x
      ylab("Total registros") + # eje y
      ggtitle("Comuna") + # titulo del grafico
      coord_flip()

# agrupar por disenio de via
acc_diseno <- acc %>% 
              group_by(DISENO) %>% 
              summarize(total_registros = n())

# crear grafica disenio 
p5 <- ggplot(data = acc_diseno, aes(x = DISENO, y = total_registros)) +
      geom_bar(stat = "identity", 
               position = "dodge", 
               fill = "dodgerblue3", 
               color = "grey48",
               alpha = 0.8) +
      xlab("Diseño de la vía") + # eje x
      ylab("Total registros") + # eje y
      ggtitle("Tipo de vía") + # titulo del grafico
      coord_flip()

# crear multigrafico
multiplot(p4, p5, cols = 2)

## grafica tipo de accidente

# agrupar por tipo de accidente
acc_clase <- acc %>% 
              group_by(CLASE) %>% 
              summarize(total_registros = n())

# crear grafica
p6 <- ggplot(data = acc_clase, aes(x = CLASE, y = total_registros)) +
      geom_bar(stat = "identity", 
               position = "dodge", 
               fill = "dodgerblue3", 
               color = "grey48", 
               alpha =0.8) +
      xlab("Tipo de accidente") + # eje x
      ylab("Total registros") + # eje y
      ggtitle("Número de accidentes por tipo de accidente") # titulo del gráfico

# visualizar grafica
p6

## grafica gravedad

# agrupar por gravedad
acc_gravedad <- acc %>% 
              group_by(GRAVEDAD) %>% 
              summarize(total_registros = n())

# crear grafica
p7 <- ggplot(data = acc_gravedad, aes(x = GRAVEDAD, y = total_registros)) + 
      geom_bar(stat = "identity",
               position = "dodge", 
               fill = "dodgerblue3", 
               color = "grey48", 
               alpha = 0.8) +
      xlab("Gravedad del accidente") + # eje x
      ylab("Total registros")+ # eje y
      # titulo del grafico
      ggtitle("Número de accidentes por gravedad del accidente") 

# visualizar grafica
p7

## grafica tipo de accidente y gravedad

# agrupar por clase y gravedad
acc_group_clase <- acc %>% group_by(CLASE, GRAVEDAD) %>% summarize(conteo = n())

# graficar
ggplot(data = acc_group_clase, aes(x = CLASE, y = conteo, fill = GRAVEDAD)) +
  geom_bar(stat = "identity", position = "dodge") +
    coord_flip()

## grafica tipo de accidente por mes

# agrupar por clase y mes
acc_group_mes <- acc %>% group_by(MES,CLASE) %>% summarize(conteo = n()) 

# graficar
ggplot(data = acc_group_mes, aes(x = CLASE, y = conteo, fill = MES)) + 
    geom_bar(stat = "identity", position = "dodge")

## grafica por tipo de accidente por anio

# agrupar por clase y periodo
acc_group_year <- acc %>% group_by(CLASE, PERIODO) %>% summarize(conteo = n()) 

# convertir periodo a categoria
acc_group_year$PERIODO <- as.factor(acc_group_year$PERIODO)

# graficar
ggplot(data = acc_group_year, aes(x = CLASE, y = conteo, fill = PERIODO)) +
    geom_bar(stat = "identity", position = "dodge") +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))+
    scale_fill_manual(values = c("#d0d1e6",
                                 "#bdc9e1",
                                 "#74a9cf",
                                 "#2b8cbe",
                                 "#045a8d"))

## grafica por tipo de accidente por dia de la semana

# agrupar por clase y dia_nombre
acc_group_dia <- acc %>% group_by(CLASE, DIA_NOMBRE) %>% summarize(conteo = n()) 

# graficar
ggplot(data = acc_group_dia, aes(x = CLASE, y = conteo, fill = DIA_NOMBRE)) + 
    geom_bar(stat = "identity", position = "dodge") +
    scale_fill_manual(values = c("#eff3ff", 
                                 "#c6dbef", 
                                 "#9ecae1",
                                 "#6baed6", 
                                 "#4292c6", 
                                 "#2171b5", 
                                 "#084594")) +
    coord_flip()

## grafica por tipo de accidente los dias festivoos

# agrupar por clase y festivo
acc_group_festivo <- acc %>% 
                     group_by(CLASE, FESTIVO) %>% 
                     summarize(conteo = n()) 

# graficar
ggplot(data = acc_group_festivo, aes(x = CLASE, y = conteo, fill = FESTIVO)) + 
    geom_bar(stat = "identity", position = "dodge")

## grafica accidentes diarios

# agrupar
acc_fecha <- acc %>% 
             group_by(FECHA) %>% 
             summarize(total_registros = n())

# graficar
plot_ly (data = subset(acc_fecha,subset = (FECHA <= "2018-12-31")),
         x = ~FECHA,
         y = ~total_registros,
         type = "scatter",
         mode = "lines",
         line = list(width = 1, color = 'rgb(205, 12, 24)'))%>%
  layout(title = 'Registros de accidentalidad',
         xaxis = list(title = "Día"),
         yaxis = list(title = "Total registros"))

## grafica accidentes por anio

# crear columna de anio
acc_fecha$year <- format(acc_fecha$FECHA,"%Y")

# graficar
plot_ly (data = subset(acc_fecha, subset = (FECHA <= "2018-12-31")),
         x = ~FECHA,
         y = ~total_registros,
         type = "scatter",
         mode = "lines",
         split = ~year,
         line = list(width = 1))%>%
  layout(title = 'Registros de accidentalidad',
         xaxis = list(title = "Día"),
         yaxis = list(title = "Total registros"))

## promedio diario

# agregar datos por promedio
aggregate(total_registros ~ year, data = acc_fecha, FUN = mean)

## grafica del promedio diario por mes y por anio

# crear columna mes
acc_fecha$Fecha <- as.Date(acc_fecha$FECHA, "%d/%m/%Y")
acc_fecha$mes <- format(acc_fecha$Fecha, "%m")
acc_fecha$mes <- strftime(acc_fecha$Fecha, format = "%B")

# ordenar columna mes 
acc_fecha$mes <- ordered(acc_fecha$mes, levels = c( "enero", "febrero", "marzo",
                                                    "abril", "mayo","junio",
                                                    "julio", "agosto",
                                                    "septiembre", "octubre",
                                                    "noviembre", "diciembre"))

# agregar datos por promedio mes y anio y graficar
aggregate(total_registros ~ year*mes, data = acc_fecha, FUN = mean) %>%
  plot_ly(x = ~mes,
          y = ~total_registros,
          type = "scatter",
          mode = "lines",
          split = ~year,
          line = list(width = 1)) %>%
  layout(title = 'Promedio diario mensual de accidentes registrados',
         xaxis = list(title="Día"),
         yaxis = list(title="Total registros"))

## grafica de caja y bigotes por anio

# graficar
plot_ly (data = subset(acc_fecha, subset = (Fecha <= "2018-12-31")),
         x = ~year,
         y = ~total_registros,
         type = "box")%>%
  layout(title = 'Registros de accidentes',
         xaxis = list(title = "Año"),
         yaxis = list(title = "Total registros"))

## grafica de caja y bigotes por mes

# crear columna de mes
acc_fecha$diames <- format(acc_fecha$Fecha, "%d")

# graficar
plot_ly (data = subset(acc_fecha, subset = (Fecha <= "2018-12-31")),
         x = ~mes,
         y = ~total_registros,
         type = "box")%>%
  layout(title = 'Registros de accidentes',
         xaxis = list(title = "Mes"),
         yaxis = list(title = "Total registros"))

## grafica de caja y bigotes por dia de semana

# crear dia de semana
acc_fecha$dia_semana <- weekdays(acc_fecha$Fecha)
acc_fecha$dia_semana <- ordered(acc_fecha$dia_semana, 
                                levels = c( "lunes", "martes", "miércoles",
                                            "jueves", "viernes", "sábado",
                                            "domingo"))

# graficar
plot_ly (data = subset(acc_fecha, subset = (Fecha <= "2018-12-31")),
         x = ~dia_semana,
         y = ~total_registros,
         type = "box")%>%
  layout(title = 'Registros de accidentes',
         xaxis = list(title = "Mes"),
         yaxis = list(title = "Total registros"))

2.3.3. Datos atípicos o nulos

El manejo de datos nulos se dio en la sección 2.3.1. donde se tomaron las siguientes decisiones:

• Si CLASE era nulo, se eliminaba por la naturaleza del objetivo.
• Si DISENO era nulo, se eliminaba por la naturaleza del objetivo.
• Si COMUNA_BARRIO era nulo o no estaba en la tabla de barrios, se eliminaba ya que posiblemente corresponda a un dato difícil de medir.

Después de realizar el análisis, se concluye que ningún dato atípico será eliminado ya que corresponden a los datos de las tablas originales, a los cuales se les aplicaron transformaciones mínimas.

Regresar

Elección del modelo de regresión

Modelo lineal grupo 0

## regresion lineal grupo 0

# crear modelo
#acc_year$PERIODO <- as.numeric(acc_year$PERIODO)
modelo_lm_0 <- lm(ACCIDENTES ~ FECHA + CLASE + DISENO + DIA_NOMBRE + PERIODO +
                    FESTIVO + MADRE + NAVIDAD + BRUJITOS + SEMSANTA + ESCOLAR,
                  data = acc_agrupado_0, 
                  subset = (FECHA <= "2017-12-31"))

# ver resultados del modelo
summary(modelo_lm_0)

## 
## Call:
## lm(formula = ACCIDENTES ~ FECHA + CLASE + DISENO + DIA_NOMBRE + 
##     PERIODO + FESTIVO + MADRE + NAVIDAD + BRUJITOS + SEMSANTA + 
##     ESCOLAR, data = acc_agrupado_0, subset = (FECHA <= "2017-12-31"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.300  -6.657  -2.583   4.472  63.013 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          1.040e+03  5.886e+02   1.768 0.077120 .  
## FECHA                2.176e-03  7.966e-04   2.732 0.006304 ** 
## CLASECAIDA OCUPANTE -3.321e-01  2.841e-01  -1.169 0.242387    
## CLASECHOQUE          2.014e+01  2.508e-01  80.296  < 2e-16 ***
## CLASEINCENDIO       -9.749e+00  2.607e+00  -3.739 0.000185 ***
## CLASEOTRO            1.147e+00  2.814e-01   4.077 4.57e-05 ***
## CLASEVOLCAMIENTO    -5.297e+00  3.086e-01 -17.165  < 2e-16 ***
## DISENOGLORIETA      -4.280e+00  6.008e-01  -7.125 1.08e-12 ***
## DISENOINTERSECCION   9.001e+00  5.622e-01  16.011  < 2e-16 ***
## DISENOLOTE O PREDIO  6.244e+00  5.647e-01  11.057  < 2e-16 ***
## DISENOPASO A NIVEL  -5.699e+00  1.701e+00  -3.351 0.000808 ***
## DISENOPASO ELEVADO  -5.339e+00  7.015e-01  -7.611 2.86e-14 ***
## DISENOPASO INFERIOR -3.482e+00  8.257e-01  -4.217 2.48e-05 ***
## DISENOPONTON        -3.752e+00  2.581e+00  -1.454 0.146002    
## DISENOPUENTE        -5.320e+00  7.470e-01  -7.121 1.11e-12 ***
## DISENOTRAMO VIA      2.445e+01  5.506e-01  44.412  < 2e-16 ***
## DISENOTUNEL         -2.347e+00  2.653e+00  -0.885 0.376368    
## DISENOVIA PEATONAL   9.371e-01  1.756e+00   0.534 0.593660    
## DIA_NOMBREJUEVES     3.904e+00  3.139e-01  12.436  < 2e-16 ***
## DIA_NOMBRELUNES      3.618e+00  3.202e-01  11.299  < 2e-16 ***
## DIA_NOMBREMARTES     4.262e+00  3.143e-01  13.562  < 2e-16 ***
## DIA_NOMBREMIERCOLES  3.865e+00  3.152e-01  12.261  < 2e-16 ***
## DIA_NOMBRESABADO     3.258e+00  3.172e-01  10.271  < 2e-16 ***
## DIA_NOMBREVIERNES    4.386e+00  3.127e-01  14.027  < 2e-16 ***
## PERIODO             -5.409e-01  2.985e-01  -1.812 0.069971 .  
## FESTIVOSi           -2.828e+00  5.185e-01  -5.454 4.98e-08 ***
## MADRESi              1.414e+00  1.650e+00   0.857 0.391666    
## NAVIDADSi           -2.479e+00  8.407e-01  -2.949 0.003193 ** 
## BRUJITOSSi           5.316e-01  1.522e+00   0.349 0.726828    
## SEMSANTASi          -2.466e+00  6.244e-01  -3.949 7.87e-05 ***
## ESCOLARSi           -3.418e-01  6.907e-01  -0.495 0.620711    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.71 on 17203 degrees of freedom
## Multiple R-squared:  0.5463, Adjusted R-squared:  0.5455 
## F-statistic: 690.5 on 30 and 17203 DF,  p-value: < 2.2e-16

Modelo lineal grupo 1

## regresion lineal grupo 1

# crear modelo
#acc_year$PERIODO <- as.numeric(acc_year$PERIODO)
modelo_lm_1 <- lm(ACCIDENTES ~ FECHA + CLASE + DIA_NOMBRE + DIA + PERIODO +
                    FESTIVO + MADRE + NAVIDAD + BRUJITOS + SEMSANTA + ESCOLAR,
                  data = acc_agrupado_1, 
                  subset = (FECHA <= "2017-12-31"))

# ver resultados del modelo
summary(modelo_lm_1)

## 
## Call:
## lm(formula = ACCIDENTES ~ FECHA + CLASE + DIA_NOMBRE + DIA + 
##     PERIODO + FESTIVO + MADRE + NAVIDAD + BRUJITOS + SEMSANTA + 
##     ESCOLAR, data = acc_agrupado_1, subset = (FECHA <= "2017-12-31"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -60.136  -3.780  -0.428   4.098  52.536 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          2.393e+03  7.909e+02   3.026 0.002485 ** 
## FECHA                3.889e-03  1.064e-03   3.654 0.000260 ***
## CLASECAIDA OCUPANTE -1.167e+00  3.450e-01  -3.382 0.000722 ***
## CLASECHOQUE          6.563e+01  3.450e-01 190.246  < 2e-16 ***
## CLASEINCENDIO       -1.028e+01  2.277e+00  -4.516 6.41e-06 ***
## CLASEOTRO            1.156e+00  3.450e-01   3.352 0.000806 ***
## CLASEVOLCAMIENTO    -7.471e+00  3.496e-01 -21.370  < 2e-16 ***
## DIA_NOMBREJUEVES     9.262e+00  4.130e-01  22.428  < 2e-16 ***
## DIA_NOMBRELUNES      8.485e+00  4.202e-01  20.191  < 2e-16 ***
## DIA_NOMBREMARTES     9.893e+00  4.144e-01  23.870  < 2e-16 ***
## DIA_NOMBREMIERCOLES  9.359e+00  4.127e-01  22.676  < 2e-16 ***
## DIA_NOMBRESABADO     7.739e+00  4.124e-01  18.768  < 2e-16 ***
## DIA_NOMBREVIERNES    1.038e+01  4.132e-01  25.117  < 2e-16 ***
## DIA                 -1.367e-02  1.266e-02  -1.080 0.279962    
## PERIODO             -1.218e+00  4.010e-01  -3.037 0.002396 ** 
## FESTIVOSi           -6.729e+00  6.740e-01  -9.984  < 2e-16 ***
## MADRESi              4.043e+00  2.107e+00   1.919 0.055026 .  
## NAVIDADSi           -7.228e+00  1.069e+00  -6.765 1.44e-11 ***
## BRUJITOSSi           1.673e+00  2.102e+00   0.796 0.426076    
## SEMSANTASi          -6.326e+00  8.072e-01  -7.837 5.28e-15 ***
## ESCOLARSi           -1.805e+00  9.617e-01  -1.877 0.060561 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.323 on 7225 degrees of freedom
## Multiple R-squared:  0.897,  Adjusted R-squared:  0.8967 
## F-statistic:  3145 on 20 and 7225 DF,  p-value: < 2.2e-16

Al verificar los valores P, se opta por retirar la variable BRUJITOS.

## regresion lineal grupo 1 sin brujitos

# crear modelo
#acc_year$PERIODO <- as.numeric(acc_year$PERIODO)
modelo_lm_1 <- lm(ACCIDENTES ~ FECHA + CLASE + DIA_NOMBRE + DIA + PERIODO +
                    FESTIVO + MADRE + NAVIDAD + SEMSANTA + ESCOLAR,
                  data = acc_agrupado_1, 
                  subset = (FECHA <= "2017-12-31"))

# ver resultados del modelo
summary(modelo_lm_1)

## 
## Call:
## lm(formula = ACCIDENTES ~ FECHA + CLASE + DIA_NOMBRE + DIA + 
##     PERIODO + FESTIVO + MADRE + NAVIDAD + SEMSANTA + ESCOLAR, 
##     data = acc_agrupado_1, subset = (FECHA <= "2017-12-31"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -60.147  -3.778  -0.431   4.090  52.536 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          2.426e+03  7.898e+02   3.071 0.002141 ** 
## FECHA                3.934e-03  1.063e-03   3.701 0.000216 ***
## CLASECAIDA OCUPANTE -1.167e+00  3.450e-01  -3.382 0.000722 ***
## CLASECHOQUE          6.563e+01  3.450e-01 190.251  < 2e-16 ***
## CLASEINCENDIO       -1.029e+01  2.277e+00  -4.517 6.36e-06 ***
## CLASEOTRO            1.156e+00  3.450e-01   3.352 0.000806 ***
## CLASEVOLCAMIENTO    -7.471e+00  3.496e-01 -21.370  < 2e-16 ***
## DIA_NOMBREJUEVES     9.262e+00  4.129e-01  22.429  < 2e-16 ***
## DIA_NOMBRELUNES      8.493e+00  4.201e-01  20.218  < 2e-16 ***
## DIA_NOMBREMARTES     9.901e+00  4.143e-01  23.897  < 2e-16 ***
## DIA_NOMBREMIERCOLES  9.359e+00  4.127e-01  22.677  < 2e-16 ***
## DIA_NOMBRESABADO     7.747e+00  4.122e-01  18.794  < 2e-16 ***
## DIA_NOMBREVIERNES    1.039e+01  4.130e-01  25.146  < 2e-16 ***
## DIA                 -1.282e-02  1.261e-02  -1.016 0.309497    
## PERIODO             -1.234e+00  4.004e-01  -3.082 0.002062 ** 
## FESTIVOSi           -6.732e+00  6.740e-01  -9.988  < 2e-16 ***
## MADRESi              4.050e+00  2.107e+00   1.922 0.054634 .  
## NAVIDADSi           -7.239e+00  1.068e+00  -6.775 1.34e-11 ***
## SEMSANTASi          -6.327e+00  8.072e-01  -7.839 5.20e-15 ***
## ESCOLARSi           -1.809e+00  9.617e-01  -1.882 0.059933 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.322 on 7226 degrees of freedom
## Multiple R-squared:  0.897,  Adjusted R-squared:  0.8967 
## F-statistic:  3310 on 19 and 7226 DF,  p-value: < 2.2e-16

Como no se observa ningún cambio en el R cuadrado ajustado, se opta por retirar DIA.

## regresion lineal grupo 1 sin brujitos ni dia

# crear modelo
#acc_year$PERIODO <- as.numeric(acc_year$PERIODO)
modelo_lm_1 <- lm(ACCIDENTES ~ FECHA + CLASE + DIA_NOMBRE + PERIODO + 
                    FESTIVO + MADRE + NAVIDAD + SEMSANTA + ESCOLAR,
                  data = acc_agrupado_1, 
                  subset = (FECHA <= "2017-12-31"))

# ver resultados del modelo
summary(modelo_lm_1)

## 
## Call:
## lm(formula = ACCIDENTES ~ FECHA + CLASE + DIA_NOMBRE + PERIODO + 
##     FESTIVO + MADRE + NAVIDAD + SEMSANTA + ESCOLAR, data = acc_agrupado_1, 
##     subset = (FECHA <= "2017-12-31"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -60.261  -3.762  -0.431   4.105  52.692 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          2.347e+03  7.860e+02   2.986 0.002838 ** 
## FECHA                3.825e-03  1.057e-03   3.617 0.000300 ***
## CLASECAIDA OCUPANTE -1.167e+00  3.450e-01  -3.382 0.000722 ***
## CLASECHOQUE          6.563e+01  3.450e-01 190.250  < 2e-16 ***
## CLASEINCENDIO       -1.028e+01  2.277e+00  -4.515 6.43e-06 ***
## CLASEOTRO            1.156e+00  3.450e-01   3.352 0.000807 ***
## CLASEVOLCAMIENTO    -7.471e+00  3.496e-01 -21.371  < 2e-16 ***
## DIA_NOMBREJUEVES     9.260e+00  4.129e-01  22.424  < 2e-16 ***
## DIA_NOMBRELUNES      8.485e+00  4.200e-01  20.203  < 2e-16 ***
## DIA_NOMBREMARTES     9.900e+00  4.143e-01  23.894  < 2e-16 ***
## DIA_NOMBREMIERCOLES  9.359e+00  4.127e-01  22.677  < 2e-16 ***
## DIA_NOMBRESABADO     7.749e+00  4.122e-01  18.798  < 2e-16 ***
## DIA_NOMBREVIERNES    1.038e+01  4.130e-01  25.143  < 2e-16 ***
## PERIODO             -1.194e+00  3.985e-01  -2.997 0.002735 ** 
## FESTIVOSi           -6.684e+00  6.723e-01  -9.941  < 2e-16 ***
## MADRESi              4.108e+00  2.106e+00   1.951 0.051131 .  
## NAVIDADSi           -7.283e+00  1.068e+00  -6.822 9.68e-12 ***
## SEMSANTASi          -6.341e+00  8.071e-01  -7.857 4.49e-15 ***
## ESCOLARSi           -1.715e+00  9.572e-01  -1.792 0.073162 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.322 on 7227 degrees of freedom
## Multiple R-squared:  0.8969, Adjusted R-squared:  0.8967 
## F-statistic:  3494 on 18 and 7227 DF,  p-value: < 2.2e-16

Como no se observa ningún cambio en el R cuadrado ajustado, se opta por retirar MADRE.

## regresion lineal grupo 1 sin brujitos, dia, madre

# crear modelo
#acc_year$PERIODO <- as.numeric(acc_year$PERIODO)
modelo_lm_1 <- lm(ACCIDENTES ~ FECHA + CLASE + DIA_NOMBRE + PERIODO + FESTIVO +
                    NAVIDAD + SEMSANTA + ESCOLAR,
                  data = acc_agrupado_1, 
                  subset = (FECHA <= "2017-12-31"))

# ver resultados del modelo
summary(modelo_lm_1)

## 
## Call:
## lm(formula = ACCIDENTES ~ FECHA + CLASE + DIA_NOMBRE + PERIODO + 
##     FESTIVO + NAVIDAD + SEMSANTA + ESCOLAR, data = acc_agrupado_1, 
##     subset = (FECHA <= "2017-12-31"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -60.249  -3.763  -0.433   4.122  52.693 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          2.308e+03  7.859e+02   2.937 0.003328 ** 
## FECHA                3.770e-03  1.057e-03   3.566 0.000365 ***
## CLASECAIDA OCUPANTE -1.167e+00  3.451e-01  -3.382 0.000725 ***
## CLASECHOQUE          6.563e+01  3.450e-01 190.213  < 2e-16 ***
## CLASEINCENDIO       -1.029e+01  2.278e+00  -4.518 6.34e-06 ***
## CLASEOTRO            1.156e+00  3.450e-01   3.351 0.000809 ***
## CLASEVOLCAMIENTO    -7.470e+00  3.496e-01 -21.364  < 2e-16 ***
## DIA_NOMBREJUEVES     9.180e+00  4.110e-01  22.337  < 2e-16 ***
## DIA_NOMBRELUNES      8.406e+00  4.181e-01  20.105  < 2e-16 ***
## DIA_NOMBREMARTES     9.819e+00  4.123e-01  23.814  < 2e-16 ***
## DIA_NOMBREMIERCOLES  9.279e+00  4.108e-01  22.590  < 2e-16 ***
## DIA_NOMBRESABADO     7.669e+00  4.103e-01  18.693  < 2e-16 ***
## DIA_NOMBREVIERNES    1.030e+01  4.110e-01  25.069  < 2e-16 ***
## PERIODO             -1.175e+00  3.985e-01  -2.948 0.003212 ** 
## FESTIVOSi           -6.691e+00  6.724e-01  -9.951  < 2e-16 ***
## NAVIDADSi           -7.299e+00  1.068e+00  -6.836 8.82e-12 ***
## SEMSANTASi          -6.358e+00  8.072e-01  -7.877 3.83e-15 ***
## ESCOLARSi           -1.711e+00  9.574e-01  -1.787 0.074014 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.324 on 7228 degrees of freedom
## Multiple R-squared:  0.8969, Adjusted R-squared:  0.8966 
## F-statistic:  3698 on 17 and 7228 DF,  p-value: < 2.2e-16

Como no se observa ningún cambio en el R cuadrado ajustado, se opta por retirar ESCOLAR.

## regresion lineal grupo 1 sin brujitos, dia, madre, escolar

# crear modelo
modelo_lm_1 <- lm(ACCIDENTES ~ FECHA + CLASE + DIA_NOMBRE + PERIODO + FESTIVO + 
                    NAVIDAD + SEMSANTA,
                  data = acc_agrupado_1, 
                  subset = (FECHA <= "2017-12-31"))

# ver resultados del modelo
summary(modelo_lm_1)

## 
## Call:
## lm(formula = ACCIDENTES ~ FECHA + CLASE + DIA_NOMBRE + PERIODO + 
##     FESTIVO + NAVIDAD + SEMSANTA, data = acc_agrupado_1, subset = (FECHA <= 
##     "2017-12-31"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -60.242  -3.769  -0.465   4.134  52.735 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          2.158e+03  7.815e+02   2.761 0.005778 ** 
## FECHA                3.561e-03  1.051e-03   3.389 0.000706 ***
## CLASECAIDA OCUPANTE -1.167e+00  3.452e-01  -3.381 0.000726 ***
## CLASECHOQUE          6.563e+01  3.451e-01 190.185  < 2e-16 ***
## CLASEINCENDIO       -1.027e+01  2.278e+00  -4.507 6.67e-06 ***
## CLASEOTRO            1.156e+00  3.451e-01   3.351 0.000810 ***
## CLASEVOLCAMIENTO    -7.468e+00  3.497e-01 -21.356  < 2e-16 ***
## DIA_NOMBREJUEVES     9.150e+00  4.107e-01  22.279  < 2e-16 ***
## DIA_NOMBRELUNES      8.370e+00  4.177e-01  20.039  < 2e-16 ***
## DIA_NOMBREMARTES     9.787e+00  4.120e-01  23.755  < 2e-16 ***
## DIA_NOMBREMIERCOLES  9.247e+00  4.104e-01  22.530  < 2e-16 ***
## DIA_NOMBRESABADO     7.670e+00  4.103e-01  18.693  < 2e-16 ***
## DIA_NOMBREVIERNES    1.027e+01  4.107e-01  25.010  < 2e-16 ***
## PERIODO             -1.098e+00  3.962e-01  -2.772 0.005588 ** 
## FESTIVOSi           -6.660e+00  6.723e-01  -9.906  < 2e-16 ***
## NAVIDADSi           -7.259e+00  1.068e+00  -6.799 1.13e-11 ***
## SEMSANTASi          -6.352e+00  8.073e-01  -7.869 4.11e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.326 on 7229 degrees of freedom
## Multiple R-squared:  0.8968, Adjusted R-squared:  0.8966 
## F-statistic:  3928 on 16 and 7229 DF,  p-value: < 2.2e-16

Dado que al retirar la variable anterior el R cuadrado ajustado bajó 0.0001, por lo cual se dejará el modelo con la variable ESCOLAR.

## regresion lineal grupo 1

# crear modelo
#acc_year$PERIODO <- as.numeric(acc_year$PERIODO)
modelo_lm_1 <- lm(ACCIDENTES ~ FECHA + CLASE + DIA_NOMBRE + PERIODO + FESTIVO + 
                    NAVIDAD + SEMSANTA + ESCOLAR,
                  data = acc_agrupado_1, 
                  subset = (FECHA <= "2017-12-31"))

# ver resultados del modelo
summary(modelo_lm_1)

## 
## Call:
## lm(formula = ACCIDENTES ~ FECHA + CLASE + DIA_NOMBRE + PERIODO + 
##     FESTIVO + NAVIDAD + SEMSANTA + ESCOLAR, data = acc_agrupado_1, 
##     subset = (FECHA <= "2017-12-31"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -60.249  -3.763  -0.433   4.122  52.693 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          2.308e+03  7.859e+02   2.937 0.003328 ** 
## FECHA                3.770e-03  1.057e-03   3.566 0.000365 ***
## CLASECAIDA OCUPANTE -1.167e+00  3.451e-01  -3.382 0.000725 ***
## CLASECHOQUE          6.563e+01  3.450e-01 190.213  < 2e-16 ***
## CLASEINCENDIO       -1.029e+01  2.278e+00  -4.518 6.34e-06 ***
## CLASEOTRO            1.156e+00  3.450e-01   3.351 0.000809 ***
## CLASEVOLCAMIENTO    -7.470e+00  3.496e-01 -21.364  < 2e-16 ***
## DIA_NOMBREJUEVES     9.180e+00  4.110e-01  22.337  < 2e-16 ***
## DIA_NOMBRELUNES      8.406e+00  4.181e-01  20.105  < 2e-16 ***
## DIA_NOMBREMARTES     9.819e+00  4.123e-01  23.814  < 2e-16 ***
## DIA_NOMBREMIERCOLES  9.279e+00  4.108e-01  22.590  < 2e-16 ***
## DIA_NOMBRESABADO     7.669e+00  4.103e-01  18.693  < 2e-16 ***
## DIA_NOMBREVIERNES    1.030e+01  4.110e-01  25.069  < 2e-16 ***
## PERIODO             -1.175e+00  3.985e-01  -2.948 0.003212 ** 
## FESTIVOSi           -6.691e+00  6.724e-01  -9.951  < 2e-16 ***
## NAVIDADSi           -7.299e+00  1.068e+00  -6.836 8.82e-12 ***
## SEMSANTASi          -6.358e+00  8.072e-01  -7.877 3.83e-15 ***
## ESCOLARSi           -1.711e+00  9.574e-01  -1.787 0.074014 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.324 on 7228 degrees of freedom
## Multiple R-squared:  0.8969, Adjusted R-squared:  0.8966 
## F-statistic:  3698 on 17 and 7228 DF,  p-value: < 2.2e-16

Modelo lineal grupo 2

## regresion lineal grupo 2

# crear modelo
#acc_year$PERIODO <- as.numeric(acc_year$PERIODO)
modelo_lm_2 <- lm(ACCIDENTES ~ FECHA + DIA_NOMBRE + PERIODO + CLASE + MES + 
                    DIA + COMUNA + FESTIVO + MADRE + NAVIDAD + BRUJITOS + 
                    SEMSANTA + ESCOLAR,
                  data = acc_agrupado_2, 
                  subset = (FECHA <= "2017-12-31"))

# ver resultados del modelo
summary(modelo_lm_2)

## 
## Call:
## lm(formula = ACCIDENTES ~ FECHA + DIA_NOMBRE + PERIODO + CLASE + 
##     MES + DIA + COMUNA + FESTIVO + MADRE + NAVIDAD + BRUJITOS + 
##     SEMSANTA + ESCOLAR, data = acc_agrupado_2, subset = (FECHA <= 
##     "2017-12-31"))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.906 -1.658 -0.210  1.170 35.549 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                -3.050e+04  1.252e+04  -2.436 0.014840 *  
## FECHA                      -4.232e-02  1.739e-02  -2.434 0.014936 *  
## DIA_NOMBREJUEVES            8.658e-01  4.405e-02  19.657  < 2e-16 ***
## DIA_NOMBRELUNES             8.187e-01  4.477e-02  18.285  < 2e-16 ***
## DIA_NOMBREMARTES            9.247e-01  4.417e-02  20.936  < 2e-16 ***
## DIA_NOMBREMIERCOLES         8.776e-01  4.406e-02  19.918  < 2e-16 ***
## DIA_NOMBRESABADO            7.072e-01  4.401e-02  16.068  < 2e-16 ***
## DIA_NOMBREVIERNES           1.010e+00  4.405e-02  22.936  < 2e-16 ***
## PERIODO                     1.548e+01  6.354e+00   2.436 0.014846 *  
## CLASECAIDA OCUPANTE        -1.801e-01  3.857e-02  -4.670 3.02e-06 ***
## CLASECHOQUE                 3.896e+00  3.292e-02 118.364  < 2e-16 ***
## CLASEINCENDIO              -1.104e+00  6.777e-01  -1.629 0.103225    
## CLASEOTRO                  -8.732e-02  3.742e-02  -2.333 0.019639 *  
## CLASEVOLCAMIENTO           -6.746e-01  4.952e-02 -13.623  < 2e-16 ***
## MES                         1.300e+00  5.293e-01   2.455 0.014085 *  
## DIA                         4.131e-02  1.740e-02   2.374 0.017611 *  
## COMUNAARANJUEZ              1.645e+00  1.499e-01  10.977  < 2e-16 ***
## COMUNABELEN                 2.018e+00  1.503e-01  13.430  < 2e-16 ***
## COMUNABUENOS AIRES          9.004e-01  1.514e-01   5.945 2.77e-09 ***
## COMUNACASTILLA              2.925e+00  1.490e-01  19.629  < 2e-16 ***
## COMUNADOCE DE OCTUBRE       7.325e-01  1.519e-01   4.823 1.42e-06 ***
## COMUNAEL POBLADO            2.974e+00  1.515e-01  19.640  < 2e-16 ***
## COMUNAGUAYABAL              2.303e+00  1.504e-01  15.307  < 2e-16 ***
## COMUNALA AMERICA            7.047e-01  1.529e-01   4.609 4.05e-06 ***
## COMUNALA CANDELARIA         5.855e+00  1.483e-01  39.493  < 2e-16 ***
## COMUNALAURELES              3.665e+00  1.495e-01  24.513  < 2e-16 ***
## COMUNAMANRIQUE              7.541e-01  1.517e-01   4.973 6.61e-07 ***
## COMUNAPALMITAS             -1.431e+00  9.416e-01  -1.520 0.128463    
## COMUNAPOPULAR               4.487e-01  1.570e-01   2.858 0.004264 ** 
## COMUNAROBLEDO               1.833e+00  1.495e-01  12.262  < 2e-16 ***
## COMUNASAN ANTONIO DE PRADO -3.302e-01  1.690e-01  -1.953 0.050783 .  
## COMUNASAN CRISTOBAL         1.195e-01  1.621e-01   0.737 0.461059    
## COMUNASAN JAVIER            3.406e-01  1.548e-01   2.200 0.027795 *  
## COMUNASANTA CRUZ            1.187e-01  1.582e-01   0.751 0.452747    
## COMUNASANTA ELENA           1.452e-01  1.835e-01   0.791 0.429023    
## COMUNAVILLA HERMOSA         5.328e-01  1.528e-01   3.488 0.000488 ***
## FESTIVOSi                  -6.547e-01  7.309e-02  -8.956  < 2e-16 ***
## MADRESi                     4.205e-01  2.236e-01   1.881 0.060009 .  
## NAVIDADSi                  -5.505e-01  1.145e-01  -4.806 1.54e-06 ***
## BRUJITOSSi                  6.275e-02  2.117e-01   0.296 0.766923    
## SEMSANTASi                 -5.307e-01  9.041e-02  -5.870 4.38e-09 ***
## ESCOLARSi                  -1.328e-01  1.017e-01  -1.306 0.191710    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.791 on 58670 degrees of freedom
## Multiple R-squared:  0.4389, Adjusted R-squared:  0.4385 
## F-statistic:  1119 on 41 and 58670 DF,  p-value: < 2.2e-16

Agregando la variable de comunas se observa que el modelo es menos significativo que el modelo agrupado 1, por lo que se opta por hacer la validación con el modelo agrupado 1.

Calculo de \(R^2\) manual de la regresión lineal

## calcular r2 del modelo lineal

# promedio variable observada de entrenamiento
y0_tr <- mean(acc_agrupado_1$ACCIDENTES[acc_agrupado_1$FECHA <= "2017-12-31"])

# calcular distancia variable observada al promedio 
r0_tr <- acc_agrupado_1$ACCIDENTES[acc_agrupado_1$FECHA <= "2017-12-31"] - y0_tr

# calcular la suma cuadratica total (TSS)
R0_tr<-mean(r0_tr^2)

# hallar variable estimada
y_pred_tr_lm<-predict(modelo_lm_1)

# calcular distancia variable observaa y variable estimada
r_tr_lm <- acc_agrupado_1$ACCIDENTES[acc_agrupado_1$FECHA <= "2017-12-31"] -
           y_pred_tr_lm

# calcular suma cuadratica de la regresion (ESS)
R_tr_lm<-mean(r_tr_lm^2)

# calculo del r2 (r2 = 1 - (ESS/TSS))
R2_lm<-1-R_tr_lm/R0_tr

# ver r2
print(R2_lm)

## [1] 0.8968859

Modelo Poisson

Se procede a evaluar el modelo lineal generalizado de una distribución Poisson en el grupo 1.

## regresion lineal poisson grupo 1

# crear modelo
modelo_glm <- glm(ACCIDENTES ~ FECHA + CLASE + DIA_NOMBRE + PERIODO +
                    FESTIVO + NAVIDAD + SEMSANTA + ESCOLAR,
                  data = acc_agrupado_1,
                  subset = (FECHA <= "2017-12-31"),
                  family = "poisson")

# ver resultados del modelo
summary(modelo_glm)

## 
## Call:
## glm(formula = ACCIDENTES ~ FECHA + CLASE + DIA_NOMBRE + PERIODO + 
##     FESTIVO + NAVIDAD + SEMSANTA + ESCOLAR, family = "poisson", 
##     data = acc_agrupado_1, subset = (FECHA <= "2017-12-31"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -7.6029  -1.0033  -0.0912   0.7844   6.5343  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          1.017e+02  1.752e+01   5.807 6.35e-09 ***
## FECHA                1.630e-04  2.355e-05   6.922 4.47e-12 ***
## CLASECAIDA OCUPANTE -1.084e-01  1.131e-02  -9.586  < 2e-16 ***
## CLASECHOQUE          1.916e+00  8.328e-03 230.063  < 2e-16 ***
## CLASEINCENDIO       -2.432e+00  2.427e-01 -10.020  < 2e-16 ***
## CLASEOTRO            9.722e-02  1.074e-02   9.055  < 2e-16 ***
## CLASEVOLCAMIENTO    -1.065e+00  1.564e-02 -68.081  < 2e-16 ***
## DIA_NOMBREJUEVES     4.661e-01  1.012e-02  46.048  < 2e-16 ***
## DIA_NOMBRELUNES      4.335e-01  1.035e-02  41.900  < 2e-16 ***
## DIA_NOMBREMARTES     4.898e-01  1.006e-02  48.680  < 2e-16 ***
## DIA_NOMBREMIERCOLES  4.689e-01  1.010e-02  46.447  < 2e-16 ***
## DIA_NOMBRESABADO     4.015e-01  1.021e-02  39.312  < 2e-16 ***
## DIA_NOMBREVIERNES    5.096e-01  1.002e-02  50.858  < 2e-16 ***
## PERIODO             -5.082e-02  8.881e-03  -5.722 1.05e-08 ***
## FESTIVOSi           -3.366e-01  1.736e-02 -19.392  < 2e-16 ***
## NAVIDADSi           -3.737e-01  2.858e-02 -13.074  < 2e-16 ***
## SEMSANTASi          -3.198e-01  2.094e-02 -15.275  < 2e-16 ***
## ESCOLARSi           -7.193e-02  2.094e-02  -3.435 0.000592 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 198916  on 7245  degrees of freedom
## Residual deviance:  13692  on 7228  degrees of freedom
## AIC: 45392
## 
## Number of Fisher Scoring iterations: 4

Cálculo del Pseudo \(R^2\) del modelo lineal general de Poisson

## calcular el pseudo r2 del modelo poisson

# calcular variable estimada
y_pred_tr_glm <- predict(modelo_glm, type = "response")

# calcular diferencia variable observada y variable estimada
r_tr_glm <- acc_agrupado_1$ACCIDENTES[acc_agrupado_1$FECHA <= "2017-12-31"] -
            y_pred_tr_glm

# hallar la proxy de la suma cuadratica de la regresion (ESS)
R_tr_glm <- mean(r_tr_glm^2)

# calcular el pseudo r2 (pseudor2 = 1 - (proxy(ESS)/TSS))
R2_tr_glm <- 1 - R_tr_glm/R0_tr

# ver pseudo r2
print(R2_tr_glm)

## [1] 0.9345265

Este modelo muestra un pseudo-\(R^2\) que se usará para hacer las comparaciones entre los modelos.

Análisis de predichos versus observados

Se elabora un DataFrame para poder graficar predichos y observados de los modelos.

## crear dataframe predichos y observados

# crear dataframe
resultados_lm_glm <- data.frame(
  FECHA =   acc_agrupado_1$FECHA[acc_agrupado_1$FECHA <= "2017-12-31"],
  ACCIDENTES = acc_agrupado_1$ACCIDENTES[acc_agrupado_1$FECHA<="2017-12-31"],
  pred_lm = y_pred_tr_lm,
  pred_glm = y_pred_tr_glm,
  res_lm = r_tr_lm,
  res_glm = r_tr_glm)

Gráfica de predichos y observados

## crear grafica de predichos y observados

# crear grafica
plot_ly (data = resultados_lm_glm,
         x = ~FECHA,
         y = ~ACCIDENTES,
         type = "scatter",
         mode = "lines",
         name='Real',
         line = list(width = 1, color = 'rgb(205, 12, 24)')) %>%
  add_trace(y= ~pred_lm,
            name = 'Modelo lineal general',
            line = list(width = 1, color = 'rgb(22, 96, 167)')) %>%
  add_trace(y = ~pred_glm,
            name = 'Modelo Poisson',
            line = list(width = 1, color = 'rgb(255, 51, 0)')) %>%
  layout(title = 'Registros DE ACCIDENTES ',
         xaxis = list(title = "Fecha"),
         yaxis = list(title = "ACCIDENTES"),
         legend = list(x = 0.75, y = 0.9))

Comparativo entre modelo lineal y observados

## crear grafica de predichos y observados para el modelo lineal

# crear grafica
plot_ly (data = resultados_lm_glm,
         x = ~FECHA,
         y = ~ACCIDENTES,
         type = "scatter",
         mode = "lines",
         name = 'Real',
         line = list(width = 1, color = 'rgb(205, 12, 24)')) %>%
  add_trace(y = ~pred_lm,
            name = 'Modelo lineal general',
            line = list(width = 1, color = 'rgb(22, 96, 167)')) %>%
  layout(title = 'Registros de Accidentes',
         xaxis = list(title="Fecha"),
         yaxis = list(title="Accidentes"),
         legend = list(x = 0.75, y = 0.9))

Comparativo entre modelo Poisson y observados

## crear grafica de predichos y observados para el modelo poisson

# crear grafica
plot_ly (data = resultados_lm_glm,
         x = ~FECHA,
         y = ~ACCIDENTES,
         type = "scatter",
         mode = "lines",
         name = 'Real',
         line = list(width = 1, color = 'rgb(205, 12, 24)')) %>%
  add_trace(y = ~pred_glm,
            name = 'Modelo Ajustado Poisson',
            line = list(width = 1, color = 'rgb(22, 96, 167)')) %>%
  layout(title = 'Registros de Accidentes',
         xaxis = list(title="Fecha"),
         yaxis = list(title="Accidentes"),
         legend = list(x = 0.75, y = 0.9))

Se puede observar que los modelos logran describir las observaciones sin observar ningún comportamiento extraño.

Predichos vs observados en el modelo Poisson

## crear grafica de predichos vs observados para el modelo poisson

# crear grafica
plot_ly (data = resultados_lm_glm,
         x = ~ACCIDENTES,
         y = ~pred_lm,
         text = ~FECHA,
         type = "scatter",
         mode = "markers",
         name = 'Modelo lineal general',
         marker = list(size = 3, color = 'rgb(22, 96, 167)')) %>%
  add_trace(y = ~pred_glm,
            text = ~FECHA,
            name = 'Modelo lineal Poisson',
            marker = list(size = 3, color = 'rgb(255, 51, 0)')) %>%
  add_trace(x = c(0:150), y = c(0:150),
            mode = "lines",
            text = rep(NA, 151),
            name = "Identidad") %>%
  layout(title = 'Registros ACCIDENTALIDAD',
         xaxis = list(title = "Observados"),
         yaxis = list(title = "Predichos"),
         legend = list(x = 0.75, y = 0.9))

Residuales vs observados:

## crear grafica de residuales vs observados para el modelo poisson

# crear grafica
plot_ly (data = resultados_lm_glm,
         x = ~ACCIDENTES,
         y = ~res_lm,
         type = "scatter",
         mode = "markers",
         name = 'Modelo lineal general',
         text = ~FECHA,
         marker = list(size = 3, color = 'rgb(22, 96, 167)')) %>%
  add_trace(y = ~res_glm,
            name = 'Modelo Poisson',
            text = ~FECHA,
            marker = list(size = 3, color = 'rgb(255, 51, 0)')) %>%
  layout(title = 'Registros ACCIDENTES',
         xaxis = list(title="Observados"),
         yaxis = list(title="Residuales"),
         legend = list(x = 0.75, y = 0.9))

Predicción del 2018

Se calculan las predicciones para el modelo lineal y el modelo Poisson

## calcular valores pronosticados para los datos de validacion

# crear subconjunto de los datos de validacion
datos_val <- subset(acc_agrupado_1,
                    subset = (FECHA >= "2018-01-01" & FECHA <= "2018-12-31"))

# calcular media de datos de validacion
yo_vl <- mean(datos_val$ACCIDENTES)

# calcular valores pronosticados modelo lineal
y_pred_vl_lm <- predict(modelo_lm_1, newdata = datos_val)

# calcular valores pronosticados modelo poisson
y_pred_vl_glm <- predict(modelo_glm, type = "response", newdata = datos_val)

# calcular diferencias entre observados y pronosticados en validacion
r_vl_lm <- datos_val$ACCIDENTES - y_pred_vl_lm
r_vl_glm <- datos_val$ACCIDENTES - y_pred_vl_glm

# calcular distancia entre datos observados y su media en la validacion
r0_vl <- datos_val$ACCIDENTES - yo_vl

Se calculan los \(R^2\) de predicción en la validación

## calculo del r2 en los datos de validacion

# calcular suma cuadratica de la regresion para modelo lineal y poisson
R_vl_lm <- mean(r_vl_lm^2)
R_vl_glm <- mean(r_vl_glm^2)

# calcular suma cuadratica total
R0_vl<-mean(r0_vl^2)

# calcular r2 y pseudo r2 en la validacion
R2_vl_lm <- 1-R_vl_lm/R0_vl
R2_vl_glm <- 1-R_vl_glm/R0_vl

# visualizar resultados
print(R2_vl_lm)

## [1] 0.9032644

print(R2_vl_glm)

## [1] 0.9396879

Se crea una gráfica de los residuales en validación

## grafica de los residuales de los modelos lineales y poisson 
## para los datos de validacion

# crear dataframe
resultados_lm_glm_val <- data.frame( FECHA = datos_val$FECHA,
                                     ACCIDENTES = datos_val$ACCIDENTES,
                                     pred_lm = y_pred_vl_lm,
                                     pred_glm = y_pred_vl_glm,
                                     res_lm = r_vl_lm,
                                     res_glm = r_vl_glm)

# crear grafica
plot_ly (data = resultados_lm_glm_val,
         x = ~FECHA,
         y = ~ACCIDENTES,
         type = "scatter",
         mode = "lines",
         name = 'Real',
         line = list(width = 1, color = 'rgb(205, 12, 24)')) %>%
  add_trace(y = ~pred_lm,
            name = 'Modelo lineal general',
            line = list(width = 1, color = 'rgb(22, 96, 167)')) %>%
  add_trace(y = ~pred_glm,
            name = 'Modelo Poisson',
            line = list(width = 1, color = 'rgb(255, 51, 0)')) %>%
  layout(title = 'Registros ACCIDENTES (Validación)',
         xaxis = list(title="FECHA"),
         yaxis = list(title="ACCIDENTES"),
         legend = list(x = 0.75, y = 0.9))

Modelo árbol de decisión para el grupo 1

Para entrenar árboles de decisiones se utilizará el grupo 1 y el paquete rpart.

Entrenamiento del árbol

## ajustar un modelo de arboles de decision para el grupo 1

# cargar paquete
library(rpart)

# crear el modelo
modelo_rpart <- rpart(ACCIDENTES ~ FECHA + CLASE + DIA_NOMBRE + PERIODO + 
                        FESTIVO + NAVIDAD + SEMSANTA + ESCOLAR,
                      data = acc_agrupado_1, 
                      subset = (FECHA <= "2017-12-31"))

Visualización del resultado del modelo

## visualizar arbol

# imprimir resultados
print(modelo_rpart)

## n= 7246 
## 
## node), split, n, deviance, yval
##       * denotes terminal node
## 
## 1) root 7246 6094251.0 23.119240  
##   2) CLASE=ATROPELLO,CAIDA OCUPANTE,INCENDIO,OTRO,VOLCAMIENTO 5785  155757.7  9.521175 *
##   3) CLASE=CHOQUE 1461  633246.9 76.962350  
##     6) DIA_NOMBRE=DOMINGO 209   16374.4 45.306220 *
##     7) DIA_NOMBRE=JUEVES,LUNES,MARTES,MIERCOLES,SABADO,VIERNES 1252  372468.7 82.246810 *

Se puede visualizar de una forma más cómoda para el usuario con la función rpart.plot() del paquete rpart.plot:

## visualizar arbol

# cargar paquete
library(rpart.plot)

# visualizar resultados
rpart.plot(modelo_rpart, tweak = 1.1)

Visualización de la importancia de las variables:

## visualizar el impacto de las variables en el modelo

# ver resultados del modelo
summary(modelo_rpart)

## Call:
## rpart(formula = ACCIDENTES ~ FECHA + CLASE + DIA_NOMBRE + PERIODO + 
##     FESTIVO + NAVIDAD + SEMSANTA + ESCOLAR, data = acc_agrupado_1, 
##     subset = (FECHA <= "2017-12-31"))
##   n= 7246 
## 
##           CP nsplit  rel error     xerror        xstd
## 1 0.87053297      0 1.00000000 1.00019958 0.021783726
## 2 0.04010399      1 0.12946703 0.12967587 0.004235165
## 3 0.01000000      2 0.08936304 0.08962793 0.003517612
## 
## Variable importance
##      CLASE DIA_NOMBRE 
##         96          4 
## 
## Node number 1: 7246 observations,    complexity param=0.870533
##   mean=23.11924, MSE=841.0504 
##   left son=2 (5785 obs) right son=3 (1461 obs)
##   Primary splits:
##       CLASE      splits as  LLRLLL,  improve=0.8705330000, (0 missing)
##       DIA_NOMBRE splits as  LRRRRRR, improve=0.0113070000, (0 missing)
##       FESTIVO    splits as  RL,      improve=0.0014065900, (0 missing)
##       SEMSANTA   splits as  RL,      improve=0.0008808344, (0 missing)
##       NAVIDAD    splits as  RL,      improve=0.0006652363, (0 missing)
## 
## Node number 2: 5785 observations
##   mean=9.521175, MSE=26.9244 
## 
## Node number 3: 1461 observations,    complexity param=0.04010399
##   mean=76.96235, MSE=433.4339 
##   left son=6 (209 obs) right son=7 (1252 obs)
##   Primary splits:
##       DIA_NOMBRE splits as  LRRRRRR,     improve=0.38595340, (0 missing)
##       FESTIVO    splits as  RL,          improve=0.04166331, (0 missing)
##       NAVIDAD    splits as  RL,          improve=0.03167520, (0 missing)
##       SEMSANTA   splits as  RL,          improve=0.02596220, (0 missing)
##       FECHA      < 16530.5 to the left,  improve=0.01585194, (0 missing)
## 
## Node number 6: 209 observations
##   mean=45.30622, MSE=78.34642 
## 
## Node number 7: 1252 observations
##   mean=82.24681, MSE=297.499

Se procede con el cálculo de los \(R^2\) de predicción y entrenamiento para el árbol:

## calcular el pseudo r2 del modelo de arbol para entramiento y validacion

# hallar valores estimados y predichos
y_pred_tr_rpart <- predict(modelo_rpart)
y_pred_vl_rpart <- predict(modelo_rpart,newdata = datos_val)

# calcular diferencia entre variable observada y variable estimada
r_tr_rpart <- acc_agrupado_1$ACCIDENTES[acc_agrupado_1$FECHA <= "2017-12-31"] -
              y_pred_tr_rpart
r_vl_rpart <- datos_val$ACCIDENTES - y_pred_vl_rpart

# calcular la suma cuadratica de la estimacion 
R_tr_rpart <- mean(r_tr_rpart^2)
R_vl_rpart <- mean(r_vl_rpart^2)

# calcular los pseudo r2
R2_tr_rpart <- 1 - R_tr_rpart/R0_tr
R2_vl_rpart <- 1 - R_vl_rpart/R0_vl

# ver resultados
print(R2_tr_rpart)

## [1] 0.910637

print(R2_vl_rpart)

## [1] 0.9176783

Comparación de todos los modelos

Se crea una nueva tabla que permita comparar los modelos en el entrenamiento:

## crear dataframe de resultados en entrenamiento

# crear dataframe
comparacion_tr <- data.frame(
  FECHA = acc_agrupado_1$FECHA[acc_agrupado_1$FECHA <= "2017-12-31"],
  ACCIDENTES = acc_agrupado_1$ACCIDENTES[acc_agrupado_1$FECHA <= "2017-12-31"],
  glm = y_pred_tr_glm,
  arbol = y_pred_tr_rpart)

Se grafican todas las series en entrenamiento:

## graficas de los 3 modelos en entrenamiento y datos observados

# graficar
plot_ly (data = comparacion_tr,
         x = ~FECHA,
         y = ~ACCIDENTES,
         type = "scatter",
         mode = "lines",
         name = 'Real',
         line = list(width = 1, color = 'rgb(205, 12, 24)')) %>%
  add_trace(y = ~glm,
            name = 'Modelo Poisson',
            line = list(width = 1, color = 'rgb(22, 96, 167)')) %>%
  add_trace(y = ~arbol,
            name ='Árbol',
            line = list(width = 1, color = 'rgb(255, 51, 0)')) %>%
  layout(title = 'Registros ACCIDENTES (Entrenamiento)',
         xaxis = list(title = "FECHA"),
         yaxis = list(title = "ACCIDENTES"),
         legend = list(x = 0.75, y = 0.9))

Se crea una nueva tabla que permita comparar los modelos en la validación:

## crear dataframe de resultados en validacion

# crear dataframe
comparacion_vl <- data.frame(FECHA = datos_val$FECHA,
                             ACCIDENTES = datos_val$ACCIDENTES,
                             glm = y_pred_vl_glm,
                             arbol = y_pred_vl_rpart)

Se grafican todas las series en validación:

## graficas de los 3 modelos en entrenamiento y datos observados

# graficar
plot_ly (data = comparacion_vl,
         x = ~FECHA,
         y = ~ACCIDENTES,
         type = "scatter",
         mode = "lines",
         name = 'Real',
         line = list(width = 1, color = 'rgb(205, 12, 24)')) %>%
  add_trace(y= ~glm,
            name = 'Modelo Poisson',
            line = list(width = 1, color = 'rgb(22, 96, 167)')) %>%
  add_trace(y = ~arbol,
            name = 'Árbol',
            line = list(width = 1, color = 'rgb(255, 51, 0)')) %>%
  layout(title = 'Registros ACCIDENTES (Validación)',
         xaxis = list(title = "FECHA"),
         yaxis = list(title = "ACCIDENTES"),
         legend = list(x = 0.75, y = 0.9))

Regresar